%% GLOBAL VARIABLES

% Declare global variables:

global samp
global Shocks obs_shares
global N_types l_t first_type first_wproj can_add Add_cell
global N N_Z feasible isolated residual1 residual2 Z_names
global Y_t nco_t skl_t gen_t 
global C_sparse hyp_lim

% Set LP optimization options:

global LPopts
LPopts = optimoptions('linprog', 'Algorithm','interior-point-legacy', 'MaxIterations',500, 'OptimalityTolerance',0.0001373438, 'Display','final');



%% LOAD DATA

% Load analytic file:

if samp
    load analytic_c1_samp.mat
else
    load analytic_c1_pop.mat
end

% List of individual characteristics:

Z_names = char(fieldnames(feasible));

% Change number storage:

N = double(N);
N_types = double(N_types);

for a = 1:size(Z_names,1)
	z = Z_names(a,:);
	N_Z.(z) = double( N_Z.(z) );
end

Add_1auth = int32(Add_1auth);

Add_2auth.t  = int32(Add_2auth.t );
Add_2auth.tb = int32(Add_2auth.tb);
Add_2auth.s  = int32(Add_2auth.s );

Xp_t = single(Xp_t);

N_t = single(N_t);

skl_t = single(skl_t);
gen_t = single(gen_t);

skl_t(skl_t < 0) = NaN;
gen_t(gen_t < 0) = NaN;

% Number of coauthors:

nco_t = N_t - 1;
nco_t(l_t==0, 1) = 0;

% Sparse matrix for Cholesky decomposition:

if samp
    C_sparse = sparse(C);
end

% Draw shocks, with antithetic acceleration:

rng(333);
temp = randn(30000, max(l_t));
temp = sort(temp, 2);

Shocks = [temp; -temp(:,[3 2 1])];

% Check means:

mean(Shocks)



%% ARRAYS FOR ADDING PROJECTS

% Cell array for adding projects:

Add_cell = cell(max( max(Add_1auth(:,1)), max(Add_2auth.t) ), N_types);

% Projects with two authors:

for t = unique(Add_2auth.t)'
	for tb = unique( Add_2auth.tb( Add_2auth.t == t ) )'
		Add_cell{t,tb} = unique( Add_2auth.s( Add_2auth.t == t & Add_2auth.tb == tb ) );
	end
end

% Projects with one author:
% NOTE: Combines with any existing desired types for two-author projects.

for i = 1:length(Add_1auth)
	t = Add_1auth(i,1);
	tb = Add_1auth(i,2);
	Add_cell{t,tb} = unique( cat(1, t, Add_cell{t,tb}) );
end

% Index for first non-isolated type:

first_wproj = find(l_t > 0, 1, 'first');

% Indicator for types that can add projects:
% (rows with any non-empty cells in cell array)

can_add = false(N_types,1);

temp = ~cellfun('isempty', Add_cell);

can_add(find(sum(temp,2) > 0) ) = true;



%% CHECK

% List objects:

whos

clear a i t tb z temp

% Size of indicator for adding projects:
% (MATLAB does not check that logical index is same size as array)

size(can_add) == size(types)
